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O . Abstract 

' This paper presents a computational framework for quasi-static brittle 

y^ fracture in three dimensional solids. The paper set outs the theoretical 

^ basis for determining the initiation and direction of propagating cracks 

O based on the concept of configurational mechanics, consistent with Grif- 

O fith's theory. Resolution of the propagating crack by the finite element 

C^ mesh is achieved by restricting cracks to element faces and adapting the 

^ mesh to align it with the predicted crack direction. A local mesh improve- 

C/^ ment procedure is developed to maximise mesh quality in order to improve 

^^ both accuracy and solution robustness and to remove the influence of the 

1-C initial mesh on the direction of propagating cracks. An arc-length control 

, ^"^ technique is derived to enable the dissipative load path to be traced. A 

hierarchical hp-refinement strategy is implemented in order to improve 

^— H both the approximation of displacements and crack geometry. The per- 

^ formance of this modelling approach is demonstrated on two numerical 

^^^ examples that qualitatively illustrate its ability to predict complex crack 

CO paths. All problems are three-dimensional, including a torsion problem 

^"^ that results in the accurate prediction of a doubly-curved crack. 

o 
en 



O 1 Introduction 



Fracture is a pervasive problem in materials and structures and predictive mod- 
elling of crack propagation remains one of the most significant challenges in solid 
mechanics. A computational framework for modelling crack propagation must 
not only be able to predict the initiation and the direction of cracks but also to 
^ provide a numerical setting to accurately resolve the crack path. 

In the Finite Element Method, strategies for discretization of the discontinu- 
ities can be categorized as either smeared or discrete. The former is attractive 
from the point of view that the problem can be solved within a continuum set- 
ting, without the need for approximation of discontinuities or changing mesh 
connectivity. However, as strains localize, numerical difficulties arise and reg- 
ularization is required. Discrete approaches are able to directly approximate 



macroscopic crack geometry and therefore describe fractures in a more natu- 
ral and straightforward manner in terms of displacement jumps and tractions. 
Developments include introducing embedded displacement jumps within finite 
elements via additional enhanced strain modes (for example [1]) or with en- 
richment functions in the context of the partition of unity (Poll) method, e.g. 

[2]. 

A straightforward discrete crack approach is to restrict the crack path to 
element faces. In the simplest case, the predicted path can be strongly infiuenced 
by the mesh, thereby infiuencing the crack surface area, and strongly affecting 
the amount of energy dissipated. The crack path's dependence on the mesh can 
be somewhat reduced by using very fine, unstructured meshes, but this results 
in computationally expensive analyses. It is worth noting that the authors [3, 4] 
have shown that within a cohesive crack methodology applied to heterogeneous 
microstructures, the crack propagation is largely controlled by the need for the 
mesh to resolve the heterogeneities. However, in the modelling of ideally brittle 
homogeneous materials, studied here, this is clearly not the case. 

In addition to the need for the mesh to resolve the crack, it is necessary to 
employ some rational means of determining if a crack will propagate and, if so, 
the direction of crack propagation. The approach taken in this paper is based 
on the principle of maximal energy dissipation, using configurational forces to 
determine the direction of the propagating crack front. A similar technique was 
successfully adopted by a number of authors, but here we mainly follow the 
work of Miehe and colleagues [5, 6]. 

This paper is primarily concerned with predicting crack propagation in large 
three-dimensional problems. The efficiency of such problems, with a large num- 
ber of degrees of freedom, usually requires the use of an iterative solver for 
solving the system of algebraic equations. In this case, it is important to con- 
trol element quality in order to optimise the system matrix conditioning, thereby 
increasing the computational efficiency of the solver. This can be problematic in 
methods such as PoU, where the enrichment functions not only the increases the 
bandwidth of the stiffness matrix but also the matrix conditioning deteriorates 
[24] . In this paper we will show how controlling the mesh quality improves both 
the crack path predictions and the robustness of the solution algorithm. Local 
r-adaptivity is adopted to mitigate the inffuence of the mesh. 

Two numerical examples for crack propagation are presented that demon- 
strate the ability of the formulation to accurately predict crack paths without 
bias from the underlying mesh and the inffuence of mesh adaptivity and mesh 
quality control on the solution obtained. 

2 Body and crack kinematics 

Continuous evolution of propagating cracks are analysed within the context of 
the Arbitrary Lagrangian Eulerian description. This is expressed by differen- 
tiable mappings from the reference material conffguration to both the current 
spatial conffguration and the current material conffguration. In the context of 



this paper, these mappings are utihzed to independently observe the deforma- 
tion of material in the physical space Qt and the evolution of the crack surface 
in the material space ^t, see Fig. 1. 



Spatial 
configuration 




Figure 1: Deformation and structural change of a body with a propagating 
crack. 

The material coordinates X are mapped onto the spatial coordinates x via 
the familiar deformation map cp^ so that 



(p:^t ^^t, x = (^(X,t) 



(1) 



The mapping cp must be differentiable, injective, and orientation preserving 
except at the boundary. The physical displacements are: 

u = x-X (2) 

The reference configuration describes the body before crack extension, with 
mapping H of the reference coordinates % on to the material coordinates X as: 



H:^o^=^t, X = H(x,t) (3) 

This mapping represents a material structural change, which, in the context 
of this work, is an extension of the crack due to movement of the crack front. 
We can also define a mapping $ of the reference coordinates on to the spatial 
coordinates as: 

$:^o^^., x = $(X,t) (4) 

For convenience, we define a composite space-time vector containing both a 
coordinate vector and time. For example: 



iX^t) = (Xl,X2,X3,t) 

Thus, the following derivatives can be expressed as 
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where the material and spatial displacement fields are 
W = X — % and w = x — % 
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H and h are the gradients of the material and spatial maps. Given the inverse 
of (6) as 
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It can be seen from (10) that the deformation gradient is 

Given that the physical material cannot penetrate itself or reverse the orienta- 
tion of material coordinates, we have: 



(11) 



det(F) = ^ > 
^ ' det(H) 



(12) 



In addition, from Eq. (10), we can define the velocity of a material point X and 
the time derivative of the deformation gradient of a point X = const as 



and 



u = (p = w- hH-^ W = w - FW 
F = Vxw - hH-^VxW 



(13) 
(14) 



3 Crack, crack front and boundary conditions 
in material space 

Similar to [6], the crack surface is denoted as F G ^t^ with a crack front dT. 
The two faces of the crack surface are given by r+ and F" and C is a surface 
that encircles the crack front, see Figure 2. The crack surface is obtained by 
the limits F+,F~ -^ F and \C\ -> 0. The surface F in three dimensional space 
can be parameterized by coordinates ^ and 77, such that F = T{£^,r]). Defining 
TJ"'~ and T+'~ as the tangent vectors to the surface F+'~: 



^p 



5r+'-(e,r?) 



di 



and T+'" 



ar+'-(e,r?) 



r7=const 



df] 



(15) 



^=const 







Figure 2: Crack construction. 



the crack normal vectors N+' to the surfaces F"*"' are then given as 



N+. 



X T+'- = 



Spin[T+'-]T+> 



(16) 



The Spin[-] operator is introduced for convenience in calculating derivatives. 
The crack surface area is given as: 



Ar = ^(A,t + Af) = 1 ^ ^JIN'lldCdr, 



(17) 



The change of this crack surface area in the material configuration, for contin- 
uously evolving surfaces r+ and F", is expressed as 



Ar 



2 ^ 



\Y. 






(18) 



Ar. • Wd5 



where, in the limit, Ar+,Ar- -^ Ar- This quantity describes the current 
orientation of the crack surface. The evolution of Ar is driven by physical 
considerations discussed in the following sections. It should be noted that Ar is 
a well defined vector for every point on the surface F, including the crack front 
9F, and has dimensions of m~^ . 

The surface encircling the crack front C can be parameterised by two families 
of curves, see Figure 2, C{i)t = C(^, /?)|r7=const and jC{r])n = C(^, ^)|^=const. In 
the limit |C| ^ and integrals over the crack front are given by 



lim [(-)dS:= lim / / (-)dS = [ lim / (-)dS. 



(19) 



The surface F is created by a curve >Ct(0 that is swept along the path which 
is determined by the direction of Ar- The kinematic relationship between the 
crack surface area Ar and the crack front velocity W on C is given as 



Ar = I lim / Ar • Wd5' 



Aar • WdL. 



(20) 



or 



where A^r is a kinematic state variable that defines the direction of crack 
propagation. Note that the spatial velocity field w is restricted in the usual 
manner by the Dirichlet boundary conditions and interpenetration of surfaces 
r+ and r~ is not admissible. 



4 Dissipation of energy at crack front 

We consider a thermodynamic system composed of a solid elastic body with a 
crack. Forces applied on the boundary of the elastic body d^t can do work 
on the system and all energy within the volume of the body can be used to 
do mechanical work. Exchange of energy between the crack surface and the 
exterior is restricted to the crack front. No free energy is stored at the crack 
front to do mechanical work on the elastic body. Making use of Eq. (13), the 
power of external work on the elastic body is given as: 



^ := / U'tdS = [ |w • t - W • F^tj dS 



(21) 



where t is the external traction vector. The rate of change of internal energy of 
the system can be decomposed as follows: 

^ := ^r + ^^, (22) 

where ^y is the internal crack energy and ^^^ is the internal body energy. The 
former is defined as: 

^r := lAr (23) 

where 7 is the surface energy and has dimension [Nm~'^]. Since the new crack 
surface is created by the sweep of the curve Ct along the path defined by A^r, 
the change of the crack surface internal energy (noting Eq. (20)) is expressed 
as: 

^r := :^^r, = 7^r = 7 / A^r • WdL (24) 

at Jqy 

Since we assume that there is no dissipation of energy within the volume of the 
body, the change of internal body energy can be expressed by 



^, 



^r-=^J ^{F)dV, (25) 



where ^ is the specific free energy function. Given the relation in Eq. (14) and 
that dV = Vx • WdF, Eq. (25) can be expressed as 

^^, = / {P : Vxw + E : Vx W}dy (26) 

where 

P:=^^^, 5]:=^(F)l-FTp (27) 



are the first Piola-Kirchhoff stress and Eshelby stress tensors, respectively. 
Tiierefore, tiie first law of thermodynamics, ^ = "^^r + ^^t, can be expressed 
as 

/ |w-t-W-F^t|d5= / 7Aar-WdL+/ {P : Vxw+E : VxWjdF 

(28) 
In order to get a local form of the first law, the Gauss divergence theorem 
is applied to the last integral in Eq. (28) resulting in the following expression 

/ 7Aor • WdL = / w • {Vx • P}dV + / W • {Vx • 5]}dV 

JdV J^t JSSt 

+ / w • {t - PN}d^ + / W • {F^t + ENjd^ 

- / w- lim / PNd6'+ / W- lim / ENdS'. 

(29) 

The spatial and material conservation laws of linear momentum balance, for any 
point inside the body, are expressed as: 

Vx • P = 0, Vx • 5] = (30) 

Considering admissible velocity fields and stress fields in equilibrium with ex- 
ternal forces, we obtain the following: 

/ 7Aar-WdL- / W • lim / {ENjd^ = (31) 

JdV JdV l^nKoJ^/ 

Thus, the local form of the first law is: 

W • (7Aar - Gar) = (32) 

where the configurational (or material) force at the crack front, that is driving 
the crack extension, is defined as: 



Gar = lim / ENdL (33) 

Eq. (32 ) represents the equilibrium condition for the crack front. One solution 
is trivial, i.e. that there is no crack growth (W = 0); another solution is that 
(7Aar — Gar) = 0; a third solution is that W is orthogonal to (7Aar — Gar)- 
Note that the first law only defines if the crack is in equilibrium but not how it 
will evolve. 

Since the free energy of the elastic body can be transformed into work at 
the crack front (or other forms of energy) but no energy is stored on the crack 
surface that can be used to do mechanical work on the body, the local variant 
of the second law is simply given as: 

^ := 7W • Aar > (34) 



where '2> is the dissipation of energy per unit length of the crack front and is 
equivalent to the local change of the crack surface internal energy. This equation 
expresses the constraint that a physically admissible evolution of the crack is re- 
stricted to positive crack area growth at each point of the crack front. Healing is 
not thermodynamically admissible unless some other (bio/chemo/mechanical) 
process takes place. Although the second law places restrictions on the direction 
of crack evolution, it does not determine how A^r or W evolves and therefore 
has to be supplemented by a constitutive law; this is described in the next 
section. With this constitutive law, configurational mechanics provides a ther- 
modynamically consistent framework which accounts for topological changes to 
the material body, in the form of crack propagation, and delivers the material 
displacement field W and the work conjugate configurational forces G^r- 

5 Material resistance 

A straightforward criterion for crack growth, in the spirit of Griffith, is proposed: 

(/)(Gar) = Gar • A^r - ^c/2 < (35) 

where Qc is is a material parameter specifying the critical threshold of energy 
release per unit area of the crack surface F. 

The principle of maximum dissipation for each point of the crack front states 
that, for all possible Griffith-like forces GJp that satisfy Eq. (35), for a given 
crack configuration, the dissipation ^ attains its maximum for the actual con- 
figurational force Gar- Therefore, we have 

(Gar - GJr) • W > (36) 

Using the classical method of Lagrange multipliers, we transform this into an 
unconstrained minimisation problem. The Lagrangian functional is therefore 
defined as: 

if(GSr,'i) = -G;r-W + K(/.(G;r). (37) 

and the Kuhn- Tucker conditions become 
--^ = -W+/^Aar=0, /^ > 0, (/)(GJr) < 0, and ^(/)(GJr) = 0. (38) 

In an unloading situation, i.e. ^ < 0, we get the trivial solution of no crack 
growth, 

W = 0, (39) 

where the first law (28) is automatically satisfied and the crack orientation Aar 
is determined by the current crack geometry. For loading of a point at the crack 
front, k> ^ (and ^ = 0), 

W = ^Aar, (40) 



which constrains the crack extension W to be cohinear to A^r- Referring back 
to the first law (Eq. (32)) and the crack growth criterion (Eq. (35)), we therefore 
have: 

7Aar = Gr and 27 = g^ (41) 

Applying the principle of maximum dissipation, the crack propagation direction 
is determined by the configurational force. A thermodynamically admissible 
crack propagation is only possible for a positive local change of the crack surface 
area (see Eq. (34). We note that k has dimension of length and represents the 
crack surface kinematic state variable which we can identify as 



Jdv 



At = / kdL>{) (42) 



This is in agreement with the restriction implied by the second law, Eq. (34). 
It should be noted that, for simplicity, only an isotropic crack resistance is 
considered here. However, the current methodology could be easily extended 
for anisotropic materials. 

For a general loading case, the formulation presented could result in unstable 
crack propagation. Therefore, in this work, the external force is controlled such 
that the crack surface evolves quasi-statically. 

6 Spatial and material discretization 

The finite element approximation is applied for both the material and physical 
space 

Xh = $(x)X, xh = $(x)X 

W^ = *(%)W, w^ = ^(x)w 

where superscript h indicates approximation and (~) indicates nodal values. 
Three-dimensional domains are discretised with tetrahedral finite elements with 
hierarchical basis functions of arbitrary polynomial order, following the work of 
Ainsworth and Coyle [12]. The higher-order approximations are only applied 
to displacements in the spatial configuration, whereas a linear approximation is 
used for displacements in the material space. The material and spatial gradients 
of deformation are expressed in the classical form 

Hh = Bx(x)X, hh=Bx(x)5 (44) 

and the physical deformation gradient can be expressed as: 

F^ = h^(H^)-i = Bx(X)X (45) 

6.1 Crack direction 

The normal to the discretized crack surface F^, applying the FE approximation, 
is given by (c.f. Eq (16)): 



where ^ and 77 are local coordinates of an element's triangular face on the crack 
surface. This normal is constant for a linear element and is easily calculated 
at Gauss integration points for higher-order approximations. Utilising Eq. (46) 
and with reference to Eq. (18), the approximation to the change in crack area 
can be expressed as 
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TRI JtRI 
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m\ 



Spin 



dX^ 



[ di 



Spm 

Of] 



df] 



f)*d,} 



^ LtriJtri J 



ASW 



W 



(47) 



where /\xri indicates the standard FE assembly for the triangular surfaces of 
elements. The discrete matrix Ap has dimensions of length. 

It is important to reiterate that nodes on the discretised crack surface F^ are 
allowed to move in the material configuration, but the overall body shape and 
its volume are constant in time. Furthermore, admissible changes in material 
nodal coordinates cannot change the shape of the crack faces (nor the total area 
of the crack surface), except for nodes on the crack front dV^. As a consequence, 
the matrix Ap has a number of nonzero rows equal to the number of active crack 
front nodes. The columns represents the contribution of the unit nodal material 

velocities (W) to changes of the crack surface area A^. 

Each row of matrix Ap can be interpreted as containing the components of 
a node's material velocity which lead to a change of crack area. We can define 
the vector of nodal Griffith-like forces, which represent material resistance, as 






crack, material 



(Af)' 



1 /r 
2 



gc 



(48) 



where gc is a vector of size equal to the number of nodes, with gc for each 
component. These nodal Griffith-like forces have dimensions of Newtons [N], 
which are work conjugate to the material displacement W on the crack front. 



6.2 Residuals in spatial and material space 

The residual force vector in the discretised spatial conffguration is expressed as: 



ext, space 



mt, space 



where A is an the unknown load factor and f^xt, space ^^^ ^kit space 



(49) 
are the vectors 



of external and internal forces, respectively, and deffned as follows 

4t,space:=A/ *^thd5, tt,,pace=A/ (B^f PM^ (50) 

TRI JtRI TET JTET 
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where /\xet indicates the standard FE assembly for tetrahedral elements. Sim- 
ilarly, the residual vector for nodal degrees of freedom in the material configu- 
ration is given by 

yh ._ fh _ rh /r-|\ 

material * crack, material int, material V^ / 

where f^^^ material ^^ ^^^ uodal values of configurational forces 

tt,material := A /" (Exf S^d^ (52) 

TET JTET 

We note that for a continuous elastic body comprising an homogeneous ma- 
terial, the configurational forces within the volume should be zero. For a dis- 
cretised homogeneous body, any non-zero values are an indication of a solution 
error. In [13], this is utilised to improve solution quality by displacing the 
nodal positions in the material space. However, this improvement methodology 
is numerically expensive since it introduces large non-linearities into the sys- 
tem which are difficult to tackle in three-dimensional problems and can lead to 
convergence problems. For this reason, the nodal configurational force vector 
is calculated only for nodes on the crack front and an alternative approach is 
adopted to reduce any numerical errors in the vicinity of the crack front. An hp- 
adaptivity strategy is adopted utilising the hierarchical approximation basis [12] 
described above and a simple mesh refinement technique, using edge-based sub- 
division, in the vicinity of the crack front [14]. 

Calculating the nodal configurational forces only for the nodes at the crack 
front, we can write 

(Gr)/ := (f^t,materiai)/ fo^ / G {/ : A/} is crack front node} (53) 
Finally, utilising Eq. (48), the material residual Eq. (49) is now expressed as 



■"■ material 



i(A^)^ge-G^ (54) 



6.3 Discrete crack propagation criterion 

The crack propagation criterion for a discretised problem is expressed by a scalar 
product of the material residual and the nodal Ap, resulting in 

^^ := A^rLteriai = ^A^(A^)^gc - A^G^ < 0, (55) 



which can be expressed for each crack front node / as: 



gc 



((a^)^a^)"'a^g 



< (56) 



For nodes on the crack front for which the fracture criterion is satisfied, nodes 
are doubled and selected faces split. Faces are chosen based on the direction 
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of the material force (Gp)/. A procedure for face splitting is described in the 
following section. 

Eq (56) involves a matrix inversion, where the size of the inverted matrix 
is equal to the number of nodes on the crack front. A more straightforward 
alternative is: _ _ 

This has been tested and no significant difference to Eq (56) was found. There- 
fore, Eq. (57) is used for all analyses, although a more detailed investigation is 
needed in the future. 

The criteria (56) and (57) are a generalization of that presented by Giirses 
and Miehe in [6], where, based on dimensional consistency, the magnitude of 
the configurational nodal force is divided by the average length of the edges 
adjacent to the node being split: 

<;c - ^ < 0. (58) 

This latter fracture criterion is equivalent to the one used here (Eq. (57)) only 
for nodes with two collinear edges. 

6.4 Crack Propagation and Arc-length Method 

To trace the quasi-static dissipative loading path, we assume that the crack front 
is in a state of quasi-equilibrium and satisfies the crack propagation criterion 
Eq. (57), within some tolerance: 

\,.-2mii^l<e (59) 

||(Ah),||2 

where e is a parameter that depends on the characteristic mesh size. From 
an equilibrium configuration, the crack front is extended at each node on the 
crack front by one element length in the direction of the configurational force 
Gp. This crack extension is resolved by the FE mesh using an r-adaptivity 
approach, whereby, for each node on the crack front, associated element faces 
are realigned with the predicted crack direction and a discontinuity is introduced 
by splitting the mesh at these faces. This process is repeated for all nodes on the 
crack front. The methodology for choosing an appropriate face for realignment 
and for subsequently splitting the mesh differs from that described in [6] and is 
discussed in Section 8. Inevitably, realignment of element faces will reduce the 
quality of the elements in the neighbourhood of the crack front, and this is also 
discussed in Section 8. 

For the new configuration, we adapt the external load factor A to achieve 
global equilibrium using an arc-length technique. Therefore, the system of equa- 
tions for conservation of the material and spatial momentum is supplemented 
by a load control equation in the form 

TA = 1 • Ahx„ - 1 • Ahxj+\ = 0, (60) 

12 



where all nodal contributions to the crack surface area are added with the use of 
vector 1. This approach of extending the crack and subsequently determining 
the corresponding load factor means that we can simplify the analysis technique 
by avoiding the need for developing an algorithm to deal with the non-smooth 
Kuhn- Tucker loading/unloading conditions, controlling the load such that the 
quasi-static crack propagation criterion is enforced at each step. 

6.5 Linearised system of equations 

In the proposed framework, the global equilibrium solution is obtained by the 
Newton- Raphson method (see Algorithm 1), converging when the norms of 
residuals rgpatiai, ^material and rx are less than a given tolerance. This is solved 
for material and spatial displacements, as a fully coupled problem. We apply 
a standard linearisation procedure (see for example, [13, 17] for details) to the 
material rlj^g^^^^-^^ and spatial residuals Tgp^^^^, resulting in a linear system of 
equations for iteration i and load step n + 1: 



^x-^int, space x-^ material ^x-^ext, space 

^Xint, space ^X material ^Xext, space 

dy^rx 




space 
material 



where X^+i = X^ + (^X^+\ X^+i = X^ + SSi'^^ and A^+i = A^ + SX'^^. 

7 Mesh quality control 

Extension of the crack front can result in poor quality elements in the vicinity 
of the crack front, including inverted or severely distorted elements, that can 
make the problem difficult or impossible to solve. To address this, we apply a 
similar approach to the one proposed by Scherer et al. [13]. The key challenge 
is to enforce constraints to preserve element quality for each Newton-Raphson 
iteration, without influencing the physical response. Thus, we introduce a mea- 
sure of element quality for tetrahedral elements in terms of their shape and use 
this to drive mesh improvement. Here we restrict ourselves to movement of 
nodes, with no changes to element topology. By applying this methodology we 
are able to solve large problems on parallel computers using Krylov iterative 
solvers with standard preconditioners. 

The dihedral angles formed between the faces of a tetrahedron have been 
shown to be one of the most influential properties in terms of solution accuracy 
[25]. Large dihedral angles result in interpolation errors and small dihedral 
angles affect the conditioning of the stiffness matrix. The minimum sine of the 
dihedral angles has been used in [25], however, this is a non-smooth function 
that cannot be used with a Newton-Raphson solver. To overcome this, an 
alternative volume-length measure is proposed. Although this measure does 
not directly measure dihedral angles, it has been shown to be very effective 
at eliminating poor angles, thus improving stiffness matrix conditioning and 
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Algorithm 1: Solution procedure 



Data: Model parameters and problem geometry 
initialise load step n = 0; 
while body not fully fractured do 
for all nodes on crack front do 

if crack propagation criterion (57) is violated then 
I extend crack front 

end 
end 

initialise iteration i = 0; 

update reference configuration ^q ^ ^t\ 

while not converged do 



compute tangents and residuals and 

solve linearised system of equations (see (61)); 






XI 



•^X^+1; 






1; 



set iteration i ^ i 
end 

save data at load step n; 
load step n ^ n + 1; 



end 



interpolation errors [15, 16]. As the volume-length measure is a smooth function 
of vertex positions and its gradient is straightforward and computationally cheap 
to calculate, it is ideal for the problem at hand. The volume-length quality 
measure is defined as 



g(Hh) := 6^2 



Fq det(H^) 



goKH*^), 6(Hh):= 



det(Hh) 



(62) 



where g'o? Vb and /rms,o are the element quality, element volume and root mean 
square of the element's edge lengths respectively, in the reference configuration. 
H^ is the material deformation gradient, h is measure of change of element qual- 
ity, relative to the reference configuration, and (i/rms = ^rms/^rms,o is the stretch 
of /rms,o from the reference configuration, go is normalised so that an equilat- 
eral element has quality 1 and a degenerate element (zero volume) has quality 
0. Furthermore, 6=1 corresponds to no change and 6 = is a change leading 
to a degenerate element. An element edge length in the material configuration 
is expressed as 

?,(H'^) := ./Ax7(Hh)THhAx, (63) 



14 



where A%-^ and is the distance vector of edge j in the reference configuration. 
Thus Irms is Calculated as 



\ 



1 ^ 

^ /_^ Ij = kmtBfidlrms (64) 



^■^1 



To control the quality of elements, we enforce an admissible deformation H^ 
such that 

6(H^) > 7 for 7 G (0, 1) (65) 

In practice, 0.1 < 7 < 0.5 gives good results. This constraint on b is enforced 
by applying a volume-length log-barrier function [13] defined as 

^-E^7I^-l-(^e-7) (66) 

e=0 ^ ^^ 

where B is the barrier function for the change in element quality in the material 
configuration. It can be seen that the Log-Barrier function rapidly increases 
as the quality of an element reduces, and tends to infinity when the quality 
approaches the barrier 7, thus achieving our aim of penalising the worst quality 
elements. 

In order to build a solution scheme that incorporates a stabilising force that 
controls element quality, we define a pseudo 'stress' at the element level as a 
counterpart of the first Piola-Kirchhoff stress as follows 

dB . , .^^v.jfrms,0 f b 1 



where matrix Q is defined as follows 

Q := (H*^)-^ - 1^^ AX^(Ax^r. (68) 

It is worth noting that Q should be a zero matrix for a purely volumetric change 
or rigid body movement of a tetrahedral element. 

It is now possible to compute a vector of nodal pseudo 'forces' associated 
with Q as 

4aiity = A /" aB^QclT/ (69) 



quality 

TET JtET 

where a is a solution parameter that can vary from to 1 but is typically set 
to 1. Thus, the global vector of nodal material residual forces (Eq. (51)) is now 
augmented with this additional stabilisation term fq^aiity ^^ give: 

material * crack, material int, material ^^ quality V ' ^/ 

Starting from an equilibrium state (i.e. configurational forces and griffith-like 
forces are equal and the arc- length control (Eq. (60)) is satisfied), a subsequent 
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load step will generally result in a loss of equilibrium and mesh deformation will 
occur, with changes in element quality controlled by the volume-length log- 
barrier function. Thus, for any arbitrary perturbation from equlibrium, by (5X, 
5Sl or ^A, equilibrium would be recovered by the Newton-Raphson procedure 
irrespective of a. However, if a = 0, there is no control of element quality. In 
future, it is possible for a to also reflect the magnitude of the out-of-balance 
material force, providing a damping effect where the current iteration is a long 
way from equilibrium. To date, this has not been required. 

8 Determination of critical faces and mesh re- 
alignment 




Figure 3: Crack front extension. Green element faces represent admissible crack 
extension before projection on predicted crack surface extension AF^. Blue 
arrows represent nodal configurational forces and red arrows represent nodal 
Grifhth-like forces. 

To determine which element face will be aligned with the predicted crack 
front extension and then split in order to create the crack extension, first all 
nodes on the current crack front are ordered according to the degree with which 
they violate the discrete crack propagation criterion (57). Taking each of these 
nodes in turn, we consider a vector s = Si + S2 on the crack front associated 
with node A/'c, as shown in Fig. 3. Next, the normal of the new crack extension. 
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Ar^, is determined as follows: 

N^ = Spin[G|^]s (71) 

As an example, Fig. 3 illustrates three faces associated with node Mc as well as 
the predicted crack surface AF^ defined by N^. It is important to emphasise 
that, in general, there will be more than one admissible set of faces that are 
associated with node A^c, but for simplicity only one set is highlighted in Fig. 3. 
Only edges and faces which are connected directly to node Mc are considered. 
Edges and faces which are part of the body's surface or the existing crack surface 
are excluded from these considerations as the crack can only propagate in the 
volume of the body. 

For each admissible set of faces, 'free' nodes A// can be identified that can 
be moved so that the faces are aligned with the predicted crack extension AF^. 
Each admissible set of faces is tested and that set which results in the least 
worse element quality (Eq. (62)) is aligned with AF^, nodes are duplicated and 
the faces are split. The last part of the process is mesh smoothing in the vicinity 
of the crack front, utilising the methodology described previously. 

In the authors' implementation, the above methodology is achieved by con- 
structing an undirected graph tree for all admissible crack surface extensions for 
each each node Mc- The vertices of the tree represent element edges and lines 
of the tree represent element faces. The starting vertex is edge si and the last 
vertex is edge 52 . In the case that node Mc is on the surface of the body, it is 
only associated with one edge, i.e. s = si and the last vertex of the tree can 
be any edge on the body surface connected to node Mc- Each path through the 
tree represents a possible set of faces that can be realigned and split to form the 
crack extension. 

Once the tree is built, a search over all branches is carried out, i.e. all possible 
paths which link the two adjacent crack front edges are searched for the set of 
faces that, when aligned with the predicted crack extension AF^, lead to the 
least worse quality element (using the volume-length quality measure (Eq. 62)). 
Thus, this approach is optimal in minimising any possible deterioration in mesh 
quality as a result of mesh realignment. 

Following realignment of the element faces, node Mc and the selected element 
faces are doubled, creating two new surfaces for the crack extension. Next, 
mesh smoothing is undertaken for a patch of elements in the vicinity of node 
Mc utilising the volume-length log-barrier function described in Section 7 to 
ensure mesh quality. This has been found to be computationally efficient and in 
practice its execution takes an insignificant fraction of the total solution time. 
This process was found to be very robust and works well for meshes with varying 
element size and/or with initially poor quality. It is also worth noting that 
this procedure of crack extension and mesh smoothing could be implemented 
in virtually any FE system as an isolated and autonomous procedure since it 
does not involve any enrichment nor does it change mesh connectivity and only 
involves doubling of nodes and element faces. 



17 



9 Implementation and Numerical Examples 

All problems are discretised using 3D tetrahedral elements. An adaptive mesh 
refinement strategy is adopted that includes both h-refinement, using an edge- 
based splitting algorithm [14], and p-refinement using hierarchic basis functions 
identified above [12]. First, all elements 8i adjacent to the crack front are 
selected, then a larger set of elements 82 adjacent to and including set £1 are 
selected. All elements in set £2 are subject to edge-based splitting. The spatial 
nodal positions are discretised using higher-order approximations in the vicinity 
of the crack front. The spatial nodal positions in all elements in £^2\^i are 
discretised using second-order approximation functions and elements in £1 using 
third-order. This refinement strategy is automated so that the mesh is locally 
refined while the crack is propagating. This is illustrated in Fig 4. As the crack 
front moves forward, elements are removed from set £2 and elements revert to 
their original state. 

The solution strategy presented in this paper is implemented for parallel 
shared memory computers, utilising open-source libraries. MOAB, a mesh- 
oriented database [18], is used to store mesh data, including input and output 
operations and information about mesh topology. PETSc (Portable, Extensible 
Toolkit for Scientific Computation [19]) is used for parallel matrix and vector 
operations, the solution of linear system of equations and other algebraic oper- 
ations. ParMetis [20] is used for parallel mesh partitioning, using the PETSc 
native interface. 




Figure 4: hp-refinement at crack front 

To demonstrate the proposed modelling framework, two numerical examples 
are presented: torsion test and pull-out test. For both examples, a quasi-brittle 
material is considered, i.e. concrete. The torsion test represents the analysis of 



18 




Figure 5: Geometry and relative size comparison for torsion and pull-out test 



a relatively small physical sample whereas the pull-out test represents a much 
larger problem, see Fig. 5. These two tests not only provide a means of verifying 
the numerical robustness of the presented methodology but also to illustrate the 
limits of the Griffith-like fracture theory adopted in this paper. 

9.1 Pull out test 

This problem considers the pull-out of a steel anchor embedded in a con- 
crete cylindrical block, which is surrounded by a steel ring. The problem 
was simulated by Duan & Areias, Belytschko [8, 9] and Gasser & Holzapfel 
[10]. All geometrical data has been obtained from [8, 9, 10]. Young's modulus 
E = 30000 N/mm , Poison ratio u = 0.2 and Griffith energy Gf = 0.106 N/mm. 
Following [8], the anchor itself is not explicitly modelled; instead a uniform 
prescribed vertical displacement is applied to the concrete face corresponding 
to the upper face of the anchor's plate. Only one quarter of the specimen is 
modelled. 

A coarse mesh comprising 2995 nodes, with a characteristic mesh size of 
51mm, and a fine mesh comprising 14075 nodes, with a mesh size of 29mm, 
have been used in the analyses. A neo-Hookean material was used initially, but 
both displacements and strains were found to be small and there was no notice- 
able difference compared to the solution assuming a linear elastic material was 
observed. Therefore, all results presented here are for a linear elastic material. 

Fig. 6 shows the load-displacement response for the two different discretisa- 
tions. The coarse mesh analysis was undertaken without any mesh refinement 
at the crack front, whereas the fine mesh analysis included hp-refinement as 
described earlier. It can be seen that in the latter case, the oscillations in the 
response are smaller in magnitude. 

Fig. 6, also compares the load-displacement paths for simulations of the same 
problem reported by Gasser & Holzapfel [10], Areias & Belytschko [9] and Duan 
& Belytschko [8]. All these simulations utilise a cohesive crack methodology 
together with enrichment functions to approximate displacement discontinuities. 
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Figure 6: Load-displacement path and energy dissipation for pull-out problem 



The cohesive crack methodology assumes that significant energy is dissipated 
in the fracture process zone, in the vicinity of the crack front, associated with 
meso-scale cracking and incorporated into the model by a fictitious extension of 
the crack front with cohesive forces. The work of these cohesive forces is equal to 
the work dissipated in the fracture process zone. In cohesive models, energy is 
dissipated by crack opening; in contrast, in the present work, it is assumed that 
energy is exclusively dissipated in the process of creating new crack area. The 
present approach is justified by the assumption that the volume of the fracture 
process zone (where meso-cracking takes place) is significantly smaller than the 
total volume of the body and so it can be assumed that energy dissipation is 
restricted to the crack front (note that aggregates are approximately two orders 
of magnitude smaller than the diameter of the specimen). This explains the good 
correlation of the present results with those presented by Duan & Belytschko, 
see Fig. 6. 

It is useful to observe that, in the case of Gasser & Holzapfel [10], the peak 
strength matches the results of Duan & Belytschko. In the case of Areias & 
Belytschko [9] , the results show an over prediction of strength compared to all 
other results. It is important to recall that the present model is based only 
on three material parameters, i.e. two for the elastic solid and the third is the 
Griffith energy release rate; no additional crack separation law is needed. In 
spite of this limited number of material parameters, good agreement with the 
other results is achieved. 
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(b) 



Figure 7: (a) Propagating crack with contour plot on crack surface of strain 
energy density and contour plot on vertical slice of polynomial order (blue: 
linear, white: quadratic, red: cubic), (b) Final crack surface. Follow this link 
for animation of the results for this problem. 

9.2 Torsion test 

An experimental study by Brokenshire [23] of a torsion test of a plain concrete 
notched prismatic beam (400mm x 100mm x 100mm) is considered. The exper- 
imental procedure and full details of the boundary conditions and dimensions 
are described in Jefferson et al. [22]. The notch is placed at an oblique angle 
across the beam and extends to half the depth. The beam is placed in a steel 
loading frame, supported at three corners and loaded at the fourth corner. The 
experiment used aggregates with a maximum size of 10 mm; the characteristic 
size of the specimen can be considered to be 100mm. 

The beam and steel frame are discretized using tetrahedral elements. The 
mesh consists of 29941 nodes, see Fig. 8. For this problem, three numerical 
tests with the same initial mesh are investigated. The first test was without 
adaptive mesh refinement near the crack front, the second with one level of 
mesh refinement and a quadratic approximation basis near the crack front and 
the third with two levels of h-refinement and a cubic hierarchical approximation 
basis near the crack front. In order to study the rate of convergence associated 
with refinement, the propagating crack is resolved by splitting the original mesh 
before refinement (note that in the pull-out test it was the refined mesh that 
was split). Hence, the mesh refinement improves the approximation quality of 
the displacements fields, but the resolution of the crack front geometry is the 
same for all three tests. 

Fig. 8 shows the load-displacement response for the three numerical tests. 
Good numerical convergence is observed with increasing refinement and the 
arc-length control method is able to track the dissipative loading path. Fig. 9 
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Figure 8: Load-displacement path. Animation for this test is available under 
this link. 



shows the numerically predicted geometry of the cracked specimen overlain on 
to the experimentally observed geometry, demonstrating excellent agreement. 
Despite the good qualitative predictions, the numerical analyses (Fig. 8) over 
predict the experimental ultimate load by approximately 2.5 times. This ob- 
served difference is a consequence of assuming linear elastic fracture mechanics 
for a problem where the size of the fracture process zone is significant compared 
to the size of the problem. This was not an issue in the pull-out test, since 
the characteristic size of the problem was significantly larger (see Fig. 5). The 
fracture energy is not only used in the creation of the macroscopic crack surface, 
but a significant proportion is used to create meso-cracks in the vicinity of the 
crack front (Bazant [22]). Even though the crack geometry matches the experi- 
mental results very well, see Fig. 9, the macro-crack has insufficient surface area 
to match the total amount of dissipated energy. Dissipation of energy in the 
macro-cracks is not included in the constitutive model and the meso-cracks are 
not discretised. An enhancement of the constitutive model to include cohesive 
cracking is an option for future work. 

10 Conclusions 

This paper has presented a computational model for quasi-static, three-dimensional, 
brittle fracture based on configurational mechanics. The model predicts the di- 
rection of crack extension based on the principle of energy minimization in 
conjunction with a node-based Griffith-like crack criterion. The direction of 
the crack extension is aligned with the nodal configurational forces at the crack 
front and the crack extension is resolved by locally aligning element faces and 
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Figure 9: 

[22]. 



The simulated crack surface superimposed on experimental result 



then splitting the mesh. Both h- and p-refinement are utilised in the vicinity 
of the crack front to increase the accuracy of the evolving crack front geome- 
try and the approximation of the spatial displacements. A monolithic solution 
strategy, simultaneously solving for both the material displacements (i.e. crack 
extension) and the spatial displacements, is presented. In order to track the 
dissipative loading path for the ideally brittle material, an arc-length procedure 
has been derived that is tailored for this model. 

A key element of the solution strategy is the control of mesh quality for the 
evolving problem. A pseudo 'stress' for the mesh quality measure is derived as 
a function of the material gradient of deformation. This mesh quality control is 
used locally in the vicinity of the crack front following realignment of element 
faces but it is also used globally to stabilise the Newthon-Raphson method. It 
is worth noting that this mesh smoothing method could be implemented in any 
standard FE system using the pseudo 'stress', standard shape functions and 
standard procedures for numerical integration and vector and matrix assembly. 
A method of selecting the tetrahedral element faces for crack surface approxima- 
tion is presented. This method is based on an undirected graph (tree) allowing 
for efficient and robust search among admissible sets of the element faces for 
crack extension, which leads to minimal mesh distortion. 

This formulation can be extended to include the case of anisotropic materials, 
material heterogeneity and other dissipative processes, for example to capture 
nonlinear size effect, when the characteristic size of the microstructure needs to 
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Figure 10: The simulated crack surface viewed from the side and above. 

be included in the model. To achieve this, thermodynamic constraints on the 
local state variables can be derived by applying the Coleman procedure [7] . This 
could potentially allow the model to be extended to include dissipative processes 
related to crack opening, applying the cohesive crack methodology [3, 4]. In 
such a situation, the energy dissipation associated with a macroscopic crack 
and dissipation related to meso-level cracking are quantified independently. 

The model has been implemented on parallel shared memory computers 
using libraries such as PETSc [18] and MOAB [19]. The two numerical examples 
demonstrated the model's ability to accurately simulate complex crack paths 
and to trace the dissipative load-displacement path. 
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